In [1]:
%matplotlib inline
import cclib
import numpy as np
from scipy import linalg
from matplotlib import pyplot as plt
import seaborn as sns
sns.set_style("white")

In [2]:
def modiagram(mol, frag1, frag2):
    """Return a molecular orbital diagram by projecting the molecular orbitals of a molecule onto the ones of two
    different fragments.
    
    Fragments should have the same basis, geometry and atomic ordering as the molecule.
    
    Parameters
    ----------
    mol : str
        File name for the output calculation of the molecule to be read by cclib.parser.ccopen.
    frag1 : str
        File name for the output calculation of the first fragment to be read by cclib.parser.ccopen.
    frag2 : str
        File name for the output calculation of the second fragment to be read by cclib.parser.ccopen.
    
    Returns
    -------
    mole : array_like
        Energies the molecular orbitals of the molecule.
    frag1e : array_like
        Energies the molecular orbitals of the first fragment.
    frag2e : array_like
        Energies the molecular orbitals of the second fragment.
    U : array_like
        Transformation matrix from fragment molecular orbitals to molecular orbitals of the final molecule."""
    pmol = cclib.parser.ccopen("H2.out")
    pfrag1 = cclib.parser.ccopen("H2_f1.out")
    pfrag2 = cclib.parser.ccopen("H2_f2.out")

    dmol = pmol.parse()
    dfrag1 = pfrag1.parse()
    dfrag2 = pfrag2.parse()
    
    # The variables below hold MO coefficients of the supermolecule and its fragments in terms of basis set
    # functions.
    Cm_a = dmol.mocoeffs[0]
    if len(dmol.mocoeffs) > 1:
        Cm_b = dmol.mocoeffs[1]

    Ck_a = linalg.block_diag(dfrag1.mocoeffs[0], dfrag2.mocoeffs[0])
    if len(dfrag1.mocoeffs[1]) > 1 and len(dfrag2.mocoeffs[1]) > 1:
        Ck_b = linalg.block_diag(dfrag1.mocoeffs[1], dfrag2.mocoeffs[1])

    # Sm below holds the basis set overlaps, taken from the molecule calculation. The overlaps taken from the
    # fragment calculations would be
    #
    # Sk = linalg.block_diag(dfrag1.aooverlaps, dfrag2.aooverlaps)
    #
    # but this matrix would lack many cross-terms.
    Sm = dmol.aooverlaps
    
    # The following would produce a block matrix [[0, A], [A, 0]]. This shows we have the right basis and the
    # right ordering of it.
    #
    # plt.matshow(np.abs(Sk - Sm))
    # plt.colorbar()

    # Next, a matrix will be calculated that holds the transformation from the MOs of fragments to the MOs of the
    # molecule.
    Cf_a = np.linalg.multi_dot([Cm_a, Sm, Ck_a.transpose()])
    # if Ck_b and Cm_b are defined:
    #     Cf_b = np.linalg.multi_dot([Cm_b, Sm, Ck_b.transpose()])

    # The following code depicts the matrix and shows that the resulting MO orbitals are normalized:
    #
    # plt.matshow(Cf_a)
    # plt.colorbar()
    # print(np.sum(Cf_a**2, axis=0))

    return dmol.moenergies, dfrag1.moenergies, dfrag2.moenergies, Cf_a

In [3]:
def plot_modiagram(mole, frag1e, frag2e, U, thresh=.0817, mowidth=1., mosep=.5, cog=0.):
    """Plot a molecular orbital diagram.
    
    Parameters
    ----------
    thresh : float
        Minimum value for drawing dashed lines between fragment and molecule molecular orbitals.
    mowidth : float
        Width in the x-axis of each plotted molecular orbital.
    mosep : float
        Distance in the x-axis of minimum separation between molecular orbitals.
    cog : float
        x position for the center of the graph."""
    # We check first multiplicities.
    # https://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.unique.html
    # umos = [mos[0]]
    # dmos = default{0: 1}
    # last_e = mos[0]
    # for e in mos[1:]:
    #     if e == last_e:
    #         if 

    def _plot_levels(levels, mowidth=mowidth, mosep=mosep, cog=cog):
        for e in levels:
            # Write a line at energy e with width mowidth
            plt.plot((cog - (mowidth - mosep)/2., cog + (mowidth - mosep)/2.), (e, e), 'k-')

    _plot_levels(mole, cog=cog)
    _plot_levels(frag1e, cog=cog - mowidth)
    _plot_levels(frag2e, cog=cog + mowidth)
    
    frag1n = len(frag1e)
    frag2n = len(frag2e)
    for i, e in enumerate(mole):
        for j in range(len(U[i, :frag1n])):
            csqr = U[i, j]**2
            if csqr > thresh:
                begin = np.array([cog - mowidth + (mowidth - mosep)/2., frag1e[j]])
                end = np.array([cog - (mowidth - mosep)/2., e])
                
                plt.plot((begin[0], end[0]), (begin[1], end[1]), 'k--')
                plt.annotate('{:.0f}%'.format(csqr*100.), xy=(begin + end)/2.)
                
        for j in range(len(U[i, frag1n:(frag1n + frag2n)])):
            csqr = U[i, j]**2
            if csqr > thresh:
                begin = np.array([cog + (mowidth - mosep)/2., e])
                end = np.array([cog + mowidth - (mowidth - mosep)/2., frag1e[j]])
                
                plt.plot((begin[0], end[0]), (begin[1], end[1]), 'k--')
                plt.annotate('{:.0f}%'.format(csqr*100.), xy=(begin + end)/2.)

In [4]:
ethresh = 5.

mole, frag1e, frag2e, U = modiagram("H2.out", "H2_f1.out", "H2_f2.out")
mole_a = mole[0][mole[0] < ethresh]
frag1e_a = frag1e[0][frag1e[0] < ethresh]
frag2e_a = frag2e[0][frag2e[0] < ethresh]

plot_modiagram(mole_a, frag1e_a, frag2e_a, U)

plt.ylabel("Energy (eV)")
plt.xticks([])


[ORCA H2.out WARNING] Modes corresponding to rotations/translations may be non-zero.
Out[4]:
([], <a list of 0 Text xticklabel objects>)